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A  five-year  multi-tasked  research  project  was  pursued  by  the  present 
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interested  researchers  involved  in  computational  fluid  dynamics  (CFD) 
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fully  developed  three-dimensional  flow  in  curved  ducts,  a  parabolized 
Navier 'Stokes  analysis  for  developing  flow  in  curved  ducts,  an  unsteady 
Navier-Stokes  analysis  for  internal  and  external  flows,  adaptive  grid 
generation  for  one-  and  two-dimensional  viscous  flows,  analysis  of  the 
Neumann  problem  in  generalized  orthogonal  coordinates,  efficient  semi- 
implicit  solution  techniques  consisting  of  the  alternating-direction 
implicit  multi-grid  ( ADI-MG)  and  strongly  implicit  multi-grid  (SI-MG) 
methods,  the  direct  block  Gaussian  elimination  (BGE)  method  for  solution  of 
the  Poisson  equation  in  generalized  orthogonal  coordinates  and  the  ADI -BGE 
and  SI-BGE  methods  for  the  unsteady  Navier-Stokes  analysis  of  incompressible 
flows.  For  the  flow  inside  a  shear-driven  cavity,  the  asymptotic  flow  in 
curved  ducts  and  the  flow  in  doubly- inf  ini te  backstep  channel,  the  predicted 
results  provided  clarity  for  interpretation  of  the  available  corresponding 
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problems.  The  adaptive-grid  generation  procedure  developed  is  unique  and 
effectively  treats  multiple  critical  regions  of  high  gradients  typical  of 
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SECTION  1 


OBJECTIVES 

The  development  of  computational  fluid  dynamics  (CFD)  analyses  for 
complex  viscous  interacting  flows  was  pursued  by  the  present  investigators 
under  AFOSR  sponsorship  during  March  1980  -  February  1985.  The  broad 
objective  was  to  gain  a  better  understanding  of  the  basic  fluid  flow 
phenomena  present  in  complex  multi-dimensional  internal  viscous  flows.  To 
realize  this  objective,  the  approach  used  was  to  study,  initially,  isolated 
problem  areas  and  flow  features  through  the  simulation  of  model  flow 
problems.  Subsequently,  as  the  necessary  information  was  either  generated 
by  the  present  investigators  or  made  available  through  the  efforts  of  other 
researchers,  more  complex  model  flow  problems  were  treated  to  better 
understand  some  of  the  dominant  flow  features  in  turbomachinery  type  of  flow 
fields  and,  eventually,  develop  an  analysis  to  solve  the  three-dimensional 
flow  in  rotating  blade  passages. 

The  major  research  thrust  was  directed  in  analyzing  the  dominant 
features  of  the  flow  fields  in  diffusers,  ducts,  blade  passages,  etc. 
Significant  effort  was  directed  in  carefully  analyzing 

©  Geometrical  Complexities  and  Three-Dimensionality, 

©  Secondary  Recirculating  Flows  and  Streamwise  Separation, 

©  Unsteadiness, 

whereas  some  effort  was  also  made  in  better  understanding 
©  Turbulence  and  Compressibility. 

Several  major  CFD  analyses  were  initiated  to  study  these  flow  features; 
these  are  grouped  into  six  different  categories  and  are  briefly  outlined 
here. 

i.  ©  For  a  class  of  three-dimensional  viscous  flow  problems  with  a 
predominant  flow  direction  and  no  streamwise  separation, 
asymptotic  flow  fields  were  to  be  determined  using  the  derived 


variables  vorticity  and  velocity  (u,  V).  This  totally  independent 
formulation  was  to  be  used  to  verify,  at  least  in  part,  the 
significant  amount  of  developing-f low  predictions  determined 
earlier.  Also,  a  systematic  study  of  Dean's  instability  was  to  be 
carried  out  for  curved  square  ducts. 

o  Further,  for  the  curved  duct  flow  configuration  with  low-speed 
flow  at  relatively  high-Reynolds  number,  turbulence  was  to  be 
carefully  modelled  by  the  two-equation  (k-e)  differential  model. 
The  Wall-Function  (WF)  approach  as  well  as  the  Low  Reynolds-number 
Modelling  (LRM)  approach  were  to  be  studied  carefully  to  determine 
which  of  the  two  treatments  of  the  wall-region  lead  to  more 
accurate  and  efficient  solutions.  The  (k-e)  model  itself  needed 
further  study  in  order  to  modify  its  usage  with  flow 
configurations  encountering  secondary  flows  of  the  second  kind, 
namely,  those  generated  by  the  turbulence  itself;  a  corner-flow 
conf iguration  was  to  be  used  as  the  model  flow  problem  for  this 
purpose. 


Two-dimensional  viscous  separation  encountered  in  steady  flows 
with  adverse  pressure  gradient  was  to  be  studied  using  both  the 
Navier-Stokes  (NS)  equations  as  well  as  the  approximate  NS 
equations  referred  to  as  the  semi-elliptic  set  of  equations,  using 
primitive  variables.  An  asymmetric  channel  with  a  constriction, 
as  well  as  a  thick  plate,  were  used  as  model  problems.  Both 
laminar  as  well  as  turbulent  flow  separation  were  to  be 
investigated.  An  improved  semi-elliptic  analysis  was  to  be 
developed,  with  appropriate  treatment  of  the  inflow  and  outflow 
boundary  conditions. 


The  effect  of  unsteadiness  in  the  flow  fields  was  to  be  studied 
using  an  unsteady  Navier-Stokes  analysis.  This  analysis  was  to 
use  derived  variables,  namely,  vorticity  to  and  stream  function  41, 
and  was  to  be  applied  to  both  two-dimensional  and  axisymmetric 
geometries  such  as  the  doubly- inf  ini te  backstep  channel,  channel 
with  sudden  expansion,  pipe  with  sudden  expansion,  center-body 
combustor  geometry,  pipe  orifice,  etc. 

Surface-oriented  coordinates  were  to  be  used  so  as  to  permit 
approximations  to  the  NS  equations,  if  needed.  Analytically  as 
well  as  numerically  generated  grids  were  to  be  developed, 
depending  on  the  specific  body  shape  in  question.  An  adaptive 
grid-generation  analysis  was  to  be  developed  for  treating  multiple 
critical  regions.  The  adaptive  grid  procedure  to  be  developed  was 
to  have  been  applied  to  model  configurations. 


Accuracy  and  efficiency  of  the  numerical  analysis  developed  was  to 
be  continually  improved  by  careful  examination  of  both  the 
continuous  and  discretized  problems.  Boundary  and  initial 
conditions  were  to  be  developed  so  the  mathematical  problems  are 
well  posed  and  these  conditions  were  to  be  implemented  correctly. 
Highei — order  spline  methodology  was  to  be  developed  for  improving 
the  solution  accuracy.  The  efficiency  of  the  numerical  solutions 
was  improved  by  the  use  of  semi- impl ici t  methods  with  increasing 


2 


degree  of  implicitness  such  as  alternating  direction  implicit 
(ADI)  method  and  strongly  implicit  (SI)  procedure.  For  enhancing 
solution  convergence,  the  multi-grid  (MG)  technique  was  to  be 
coupled  with  the  existing  semi-implicit  methods  for  two- 
dimensional  flows.  The  various  operators  of  the  multi-grid 
methods  were  to  be  carefully  studied  so  as  to  lead  to  an  efficient 
solution  procedure.  The  multi-grid  method  was  to  be  implemented 
in  a  three-dimensional  PNS  analysis.  The  solution  convergence  was 
also  to  be  enhanced  by  use  of  a  local  time-stepping  method. 

Direct  solvers  available  for  solving  the  Poisson  equation  in 
uniform  Cartesian  coordinates  were  to  be  extended  for  use  with 
clustered  curvilinear  orthogonal  coordinates.  Cauchy-R iemann 
solvers  were  to  be  developed  for  the  solution  of  Poisson  equations 
in  rectangular  and  generalized  curvilinear  coordinates.  The  block 
Gaussian  elimination  (BGE)  method  was  to  be  developed  for  the 
solution  of  the  Poisson  equation  in  generalized  orthogonal 
coordinates. 


SECTION  2 


DESCRIPTION  OF  SIGNIFICANT  ACCOMPLISHMENTS 
All  of  the  areas  of  research  initiated  and  the  progress  as  well  as  the 
specific  achievements  made  in  these  studies  during  the  five-year  grant 
period  are  briefly  summarized  in  the  following  subsections. 

2 . 1  Viscous  Flow  in  Curved  Ducts:  Examination  of  Asymptotic  Flow  and 
Turbulence  Models 

A  generalized  formulation  was  developed  to  study  three-dimensional 
asymptotic  viscous  flow  in  curved  ducts  using  the  streamwise  vorticity  z, 
the  streamwise  velocity  w  and  the  cross-flow  velocities  in  terms  of  a  cross¬ 
flow  stream  function  ip.  The  analysis  readily  permits  the  consideration  of 
straight  as  well  as  curved  ducts  of  rectangular  as  well  as  polar  cross 
sections,  as  shown  in  Fig.  la.  For  these  flows,  the  similarity  parameter  of 
significance  is  the  Dean  number  K  rather  than  the  Reynolds  number  Re.  From 
the  investigators'  earlier  work  on  this  problem,  it  was  felt  that,  for 
highly  curved  configurations,  the  strong  coupling  between  the  primary  and 
the  secondary  flow  should  be  honored  by  the  numerical  solution  technique 
employed.  Consequently,  simultaneous  numerical  solutions  of  the  three 
second-order  coupled  partial  differential  equations  (PDE's)  governing  the 
flow  were  obtained  using  semi-implicit  methods.  Initially,  an  ADI  method 
was  used;  subsequently,  the  strongly  implicit  (SI)  method  was  necessary  for 
proper  convergence  of  cases  with  larger  values  of  K.  For  curved  ducts  of 
square  cross  section,  the  use  of  very  fine  grids  revealed  that  Dean's 
instability  occurs  at  K=125,  as  shown  in  Fig.  lb,  when  an  additional  pair  of 


secondary  vortices  first  makes  its  appearance.  The  PNS  analysis  of  K.  Ghia 
and  Sokhey  (1977)  and,  later,  of  U.  Ghia.  K.  Ghia  and  Goyal  (19-?9)  had 


predicted  this  instability  to  occur  at  K  =  143,  whereas  Cheng  et  al.  (1976) 
had  used  the  vortici ty-veloci ty  variables  (c,w,p)  and  predicted  this 


instability  to  occur  at  K-202.  In  the  latter  investigation,  it  was  also 
predicted  that  this  additional  pair  of  secondary  vortices  disappears  for 
K  _>  520  whereas,  in  the  present  study,  the  additional  pair  of  secondary 
vortices  persists  even  for  K  =  900.  The  significance  of  this  phenomenon  is 
that  this  second  pair  of  streamwise  vortices  creates  additional  pressure 
losses.  The  results  for  the  asymptotic  flow  in  curved  ducts  of  polar  cross 
section  are  shown  in  Fig.  1c.  The  cross-flow  streamline  contours  for 
K  =  100,  200  and  300  exhibit  no  symmetry  and  the  primary  vortex  pair  shows 
that  the  upper  vortex  is  slightly  weaker  than  the  lower  vortex.  Most  of  the 
results  for  laminar  flow  obtained  earlier  using  the  PMS  analysis  withstood 
the  test  of  comparison  with  these  asymptotic  results,  except  for  the  cross- 
flow  velocities.  In  general,  even  the  quantitative  agreement  between  the 
two  approaches  is  good,  thereby  providing  reliance  in  the  PNS  marching 
analysis.  The  detailed  asymptotic  results  for  these  flows  were  given  by 
K.  Ghia,  U.  Ghia  and  Shin  (1981).  The  fine-grid  solutions  for  this  flow 
problem  as  well  as  this  model  problem  itself  have  become  benchmarks,  as 
other  investigators  are  using  this  flow  configuration  to  validate  their 
results  and  analyses  by  comparing  with  the  present  investigators'  solutions. 

The  PNS  analysis  developed  by  U.  Ghia,  K.  Ghia  and  Goyal  (1979)  for 
laminar  flow  in  curved  ducts  was  extended  to  turbulent  flow  configurations. 
Turbulence  closure  was  achieved  by  a  two-equation  (k,e)  differential  model. 

A  comparative  study  of  the  wall-region  treatment  by  the  wall-function  (WF) 
method  and  the  Low-Reynolds  number  Modeling  (LRM)  method  was  carried  out  by 
Goyal,  K.  Ghia  and  U  Ghia  (1980)  using  a  curved  circular  pipe  configuration 
with  Re  =25,000.  Here,  Re  denotes  Reynolds  number  based  on  hydraulic 


diameter.  The  mean  flow  as  well  as  turbulence  quantities  were  examined  in 
detail  but  the  study  was  inconclusive  due  to  the  following.  The  form  of  the 
law  of  the  wall  used  did  not  account  for  the  effects  of  curvature;  on  the 
other  hand,  the  LRM  method  satisfies  the  wall-boundary  conditions  more 
accurately,  but  has  the  disadvantage  of  being  limited  to  flows  with  moderate 
Reynolds  number  and  requiring  much  finer  computational  grids.  The  detailed 
flow  results,  showing  the  effect  of  the  duct  aspect  ratio  and  curvature 
ratio  on  the  overall  flow,  were  given  by  Goyal,  K.  Ghia  and  U.  Ghia  (1981). 
Flow  results  were  obtained  for  curved  circular  pipes  with  ReQ  =  25,000  and 

236,000.  For  curved  rectangular  ducts,  results  were  obtained  for 

ReD  =  0.706x1 0^.  All  of  these  results  compare  well  with  other  existing 

experimental  and  numerical  results. 

The  90°  axial  corner  flow  configuration,  representative  of  the  flow 
near  the  junction  of  a  turbomachinery  blade  with  the  hub,  was  analyzed, 
using  a  veloci ty-vortici ty  formulation,  to  study  the  effect  of  turbulence, 
compressibility  and  mass  transfer  on  the  developing  flow  along  an  axial 
corner.  The  governing  equations  for  this  quasi-three-dimensional  flow  were 
obtained  as  the  limiting  equations  derived  from  the  general  corner-layer 
equations  formulated  earlier.  These  asymptotic  equations  were  solved  using 
a  semi-implicit  second-order  accurate  marching  scheme.  The  turbulence  was 
modeled  by  using  a  Cebeci-Smith  ( 1 97 4 )  type  two-layer  algebraic  model  in 
which  isotropy  is  assumed.  The  turbulence  stresses  were  also  modeled  using 
a  modified  form  of  the  Gessner-Emery  (1976)  anisotropic  model.  The  skin- 
friction  coefficient  obtained  using  the  two  different  turbulence  models  is 
presented  in  Fig.  2a.  The  conformity  between  these  two  sets  of  results  led 
to  the  conclusion  that  the  effect  of  anisotropy  is  not  significant  in  the 


corner-flow  configuration.  Figure  2b  shows  the  variations  in  the  streamwise 
as  well  as  the  cross-flow  velocities  due  to  suction  and  injection  at  the 
wall  for  both  laminar  as  well  as  turbulent  flow.  Additional  results  for 
this  flow  configuration  were  presented  by  Mikhail  and  Ghia  (1981)  for  a 
range  of  Mach  numbers  between  0  and  2.0,  with  adiabatic  as  well  as  heat- 
transfer  boundary  conditions  at  the  corner  walls;  the  effects  of  suction  and 
injection  were  also  included. 


Steady  laminar  incompressible  separated  flow  was  studied  by  U.  Ghia, 

K.  Ghia,  Rubin  and  Khosla  (1981).  Their  semi-elliptic  analysis  was 

formulated  using  the  primitive  variables  (V,p).  The  model  problem  was  that 
of  a  doubly  infinite  channel  with  an  asymmetric  constriction.  This  approach 
was  demonstrated  to  be  very  promising  by  comparing  the  results  with  those  of 
the  Navler-Stokes  analysis  of  the  authors.  Further,  Osswald  and  Ghia  (1981) 
used  a  totally  different  formulation  using  the  derived  variables,  namely, 
the  vorticity  w  and  the  stream  function  iJj,  and  showed  that  their  results 
agree  with  those  of  U.  Ghia  et  al.  (1981).  Hence,  further  work  was  carried 
out  using  the  semi-elliptic  analysis.  However,  with  a  larger  region  and  a 
finer  grid,  the  numerical  method  experienced  convergence  difficulties.  In 
hindsight,  it  seems  that,  if  the  semi-elliptic  analysis  had  been  used  in 
conjunction  with  a  staggered  grid  and  a  multi-grid  procedure,  the  steady 
separated  flow  could  have  been  studied  more  effectively.  Also,  it  should  be 
pointed  out  that  the  present  investigators  were  later  able  to  develop  an 
Interacting  Parabolized  Navier-Stokes  (IPNS) 

analysis  which  overcomes  these  difficulties.  This  latter  analysis  was 
developed  originally  under  NASA  sponsorship  and,  hence,  the  details  are  not 
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presented  here.  Figure  3  shows  the  typical  stream  function  contours  for 
separated  flow  computed  using  the  Navier-Stokes  analysis  of  Osswald  and 
K.  Ghia  (1981);  additional  results  are  also  given  in  this  reference.  A 
second  model  problem,  namely,  that  of  flow  past  a  class  of  two-dimensional 


blunt  bodies  was  also  initiated  to  study  flow  separation.  The  Navier-Stokes 


equations  in  conformal  coordinates  and  similarity-type  variables  were 
solved,  using  the  approximate  factorization  scheme  of  Beam  and 


Warming  (1977).  The  effect  of  Reynolds  number  on  the  length  of  the 
recirculation  region  was  studied  and  satisfactory  agreement  was  obtained 
with  the  data  of  Lane  and  Loehrke  (1980),  as  shown  in  Fig.  9a.  The  detailed 
flow  results,  including  those  for  the  axisymmetric  case,  were  given  by  Ghia 


and  Abdelhalim  (1982).  The  corresponding  turbulent  separated  flow  was 
studied  using  second-order  closure  via  the  (k-e)  two-equation  model.  The 
wall  region  was  treated  using  the  LRM  method.  Figure  9b  shows  the 


s 


distribution  of  the  turbulent  kinetic  energy  obtained  from  this  analysis  for 
a  slightly  blunt-shouldered  body  and  compared  with  the  measurements  of  Ota 
and  Narita  (1978)  for  a  completely  sharp-shouldered  thick  plate.  As 


expected,  the  computed  separation  was  milder  than  the  measured  one.  When 
the  streamwise  dimension  of  the  predicted  separation  bubble  was  scaled  up  to 
match  with  the  measured  separation  bubble,  the  computed  results  within  the 
separated-f low  region  compared  well  with  the  corresponding  data,  as  shown  in 
Fig.  9b.  The  detailed  flow  results  were  given  by  Abdelhalim,  U.  Ghia  and 
K.  Ghia  (1983) . 
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The  effect  of  unsteadiness  on  the  flow  field  was  studied  using  the 
model  problem  of  flow  over  a  backstep  in  a  doubly  infinite  channel.  The 
unsteady  Navier-Stokes  equations  were  formulated  in  terms  of  vorticity  <u  and 
stream  function  i|i.  Clustered  conformal  grids  were  used  and  nearly  optimum 
grid  distribution  was  arrived  at  to  attempt  to  honor  the  multiple  length 
scales  of  the  separated-flow  problem.  K.  Ghia,  Osswald  and  U.  Ghia  (1983) 
gave  the  detailed  flow  results  and  also  provided  comparisons  with  the  data 
of  Denham  and  Patrick  (1974)  as  well  as  Armaly  and  Durst  (1980).  Figure  5a 
shows  this  comparison  for  the  reattachment  length  of  the  primary  bubble  on 
the  lower  wall.  At  Reg  -  212,  the  calculations  show,  for  the  first  time, 

the  appearance  of  a  secondary  bubble  near  the  upper  wall  and, 
simultaneously,  a  marked  deviation  from  the  data  of  Armaly  and  Durst  (1980). 
This  information  was  communicated  to  those  authors  and  the  response  received 
was  that  spanwise  variation  had  indeed  been  observed  and,  hence,  the 
experimental  data  were  really  three-dimensional.  These  results  were 
carefully  examined  by  Osswald,  K.  Ghia  and  U.  Ghia  (1983)  who  also  obtained 
predictions  of  their  own,  as  given  in  Figs.  5b-5d;  their  deviation  from  the 
data  of  Armaly  and  Durst  (1980)  represents  the  effect  of  three- 
dimensionality  in  the  flow  field.  The  computed  results  for  a  configuration 
well  within  the  transition  regime  exhibits  a  ’near-limit-cycle'  behavior,  as 
shown  in  Fig.  6.  This  is  an  example  of  self-sustaining  oscillatory  motion, 
where  disturbance  mechanisms  are  developed  without  any  external  excitation. 
This  creates  an  instability  which  subsequently  leads  to  transition.  Figures 
7a  and  7b  show  the  time-averaged  stream-function  and  vorticity  contours, 
respectively.  The  overall  mechanism  causing  the  unsteadiness  was  also 
explained  by  the  present  investigators.  The  Strouhal  number  of  the  vortex 
shedding  is  0.38  and  the  large-scale  coherent  structure  observed  has  been 


recently  reproduced  by  Patera  (1985)  in  his  direct  numerical  simulation 


using  the  spectral-element  method;  for  a  similar  configuration,  Roose  and 
Kegelman  (1985)  also  observed  this  coherent  structure  in  their  experiments 
for  the  first  time. 

The  separated  flow  through  two-dimensional  and  axisymmetric  sudden 
expansions  was  studied  by  Osswald,  K.  Ghia  and  U.  Ghia  (1984).  The  results 
of  their  analysis  agreed  well  with  the  data  of  Durst,  Melling  and 
Whitelaw  (1974)  for  the  plane  sudden  expansion  and  with  the  analysis  of 
Macagno  and  Hung  (1967)  for  the  axisymmetric  case.  The  transient  results 
were  presented  for  the  cold  flow  in  an  axisymmetric  centerbody  combustion 
chamber.  McGreehan,  K.  Ghia,  U.  Ghia  and  Osswald  (1984)  simulated  the  flow 
inside  a  pipe  orifice  to  study  this  persistently  unsteady  flow  and  analyze 
the  nature  of  vortex  shedding.  The  collective  vortex  interaction  phenomenon 
of  Ho  (1981)  was  observed  in  this  configuration. 

The  cornerstone  of  any  internal  aerodynamic  analysis  is  the  ability  to 
predict  the  lift  and  drag  on  lifting  surfaces  such  as  airfoils,  or  cascades 
of  airfoils  where  concerned  with  turbomachinery  applications.  For 
compressor  cascades,  it  is  very  vital  to  understand  the  effect  of 
unsteadiness  on  the  flow  field  involving  phenomena  such  as  rotating  stall, 
individual  blade  stall,  stall  flutter,  compressor  surge,  etc.  Common  to  all 
of  these  physical  phenomena  is  the  occurrence  of  unsteady  flow  separation. 
When  the  flow  separates  over  a  body  surface,  a  'strong'  interaction  region 
appears  locally  where  the  pressure  field  in  the  flow  is  determined  by  the 
viscous  layer  rather  than  the  inviscid  flow.  Accurate  and  efficient  strong 
viscous-inviscid  analyses  are  still  under  development.  Therefore,  it  was 
decided  to  modify  the  objectives  for  the  last  year  of  the  proposal  and 


include  in  it  a  study  of  flow  over  an  isolated  airfoil.  At  low  speed,  the 
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local  'strong'  interaction  regions  arise  in  this  latter  flow  due  to 
boundary-layer  separation  and  rapid  flow  acceleration  immediately  aft  of  the 
trailing  edge.  The  Joukowski  airfoil,  with  its  sharp  trailing  edge,  is  an 
ideal  candidate  for  this  study,  since  it  also  permits  the  use  of  analytical 
transformation  metrics,  thereby  avoiding  the  error  incurred  due  to  the  use 
of  numerically  computed  metrics.  Therefore,  in  order  to  carefully  analyze 
persistently  unsteady  flow  with  massive  separation,  K.  Ghia,  Osswald  and  U. 
Ghia  (1985a)  extended  their  earlier  analysis  for  internal  flows  and  applied 
it  to  an  isolated  airfoil  at  high  angle  of  attack.  For  moderate  Re,  these 
authors  were  successful  in  simulating  the  flow  structure  in  the  highly 
unsteady  vortex-dominated  wake  region.  The  results  of  a  limit  cycle 
analysis  for  flow  past  a  Joukowski  airfoil  are  presented  in  Fig.  8.  Such 
limit-cycle  solutions  represent  realizations  of  a  strange  attractor,  as  all 
phase-space  trajectories  are  ultimately  attracted  to  the  time-asymptotic 
limit-cycle  solution.  The  detailed  evolution  of  the  pairwise  shedding  of 
vortices  from  the  airfoil  surface  was  discussed  by  K.  Ghia,  Osswald  and 
U.  Ghia  (1985b).  The  Strouhal  number  of  this  shedding  motion  was  S  =  0.18 
which  agrees  well  with  the  universal  wake-based  number  of  Roshko  (1954). 

2.4  Numerical  Grid  Generation 

The  problem  of  the  resolution  of  high-gradient  and  high-curvature 
regions  in  computational  fluid  dynamics  is  important  not  only  from 
considerations  of  truncation  errors  but  sometimes  simply  for  correct 
prediction  of  the  flow  details  in  these  regions.  Most  nonlinear  phenomena 
have  a  tendency  to  occur  in  very  thin  regions  which  may  or  may  not  be 
associated  with  boundaries,  as  shown  by  K.  Ghia,  U.  Ghia  and  Shin  (1981)  and 
U.  Ghia,  K.  Ghia  and  Shin  (1982).  In  the  simulation  of  viscous  flows  at 


high  Reynolds  number  in  these  thin  regions,  the  nonlinear  convective  terms 
became  large  in  comparison  with  the  viscous  terms.  In  this  circumstance,  it 
became  difficult  to  obtain  a  wiggle-free  solution  using  a  central-difference 
discretization  for  the  convective  terms.  A  new  analysis  was  developed  for 
generating  an  adaptive  grid  which  took  into  account  the  influence  of  the 
problem  geometry,  the  prevailing  physical  flow  phenomena,  the  flow 
parameters,  as  well  as  the  discretization  parameters  of  the  problem.  The 
grid  constraint  used  minimized  the  coefficient  of  the  convective  terms  in 
the  transformed  flow  equations  in  a  rational  manner  and  led  to  grid 
equations  which  were  analogous  to  the  inverted  form  of  the  Poisson  equations 
used  in  elliptic  grid-generation  techniques.  The  method  was  demonstrated 
with  the  help  of  nonlinear  one-dimensional  as  well  as  two-dimensional  model 
problems  by  comparing  the  predicted  solutions  with  the  corresponding  exact 
solutions.  A3  shown  in  Fig.  9a,  for  the  flow-dependent  adaptive  grid,  the 
desired  number  of  grid  points  have  migrated  to  the  region  of  the  high 
gradient  so  as  to  limit  the  magnitude  of  the  truncation  errors.  Figures  9b 
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and  9c  show  the  computed  and  analytical  results  for  Re  up  to  1 0  for  the 
nonlinear  viscous  Burgers'  equation  and  a  model  internal  flow  problem, 
respectively,  and  the  agreement  with  the  exact  analytical  solution  is 
excellent;  to  the  authors'  knowledge,  these  calculations  were  the  first  of 
their  kind.  The  detailed  results  for  a  few  one-dimensional  and  two- 
dimensional  model  problems  were  given  by  K.  Ghia,  U.  Ghia  and  Shin  (1983). 


2.5  Semi-Implicit  Numerical  Methods:  Accuracy  and  Efficiency 


For  improving  the  accuracy  of  a  given  numerical  method,  it  was  decided, 


early  in  this  research  program,  to  seek  solutions  with  accuracy  higher  than 


second  order  by  using  spline  methodology,  thereby  also  minimizing  the  number 


of  grid  points  needed  for  obtaining  solutions  for  viscous  flow  problems  with 
high-gradient  regions.  A  subsequent  period  of  this  research  program  saw  the 
development  of  a  multi-grid  (MG)  method  which  improves  the  solution  accuracy 
by  utilizing  very  fine  grids  while  still  retaining  rapid  convergence 
behavior.  This  MG  method  as  well  as  the  spline  method  are  discussed  next. 

It  was  shown  by  K.  Ghia,  Shin  and  (J.  Ghia  (1979),  using  primitive  variables 

(V.p),  and  by  Shin,  K.  Ghia  and  U.  Ghia  (1981),  using  derived  variables, 
that  for  low-speed  viscous  flow  problems  with  localized  high-gradient 
regions,  although  the  resulting  solutions  themselves  may  be  smooth,  their 
first-  and/or  second-order  derivatives  were  frequently  not  smooth.  This  led 
to  the  development  of  a  new  spline  technique  employing  a  quartic  spline 
polynomial  S(4,2)  of  deficiency  two,  i.e.,  with  two  continuous  derivatives. 
The  integrated  form  of  the  governing  equation,  which  is  generally  used  in 
finite-volume  techniques,  was  employed  in  this  analysis  to  complete  the 
equation  set.  Some  detailed  flow  results  obtained  using  the  spline  S(4,2) 
method  were  given  by  Turner,  K.  Ghia,  U.  Ghia  and  Keith  (1982).  The 
accuracy  of  spline  S(4,2)  is  comparable  to,  if  not  better  than,  that  of  the 
fourth-order  box  scheme  used  by  Wornom  (1977)  and  the  compact  differencing 
scheme  of  Hirsh  (1975).  In  a  review  of  higher-order  methods  for  the 
solution  of  incompressible  viscous  flows,  K.  Ghia  and  U.  Ghia  (1982)  had 
suggested  spline  3(4,2)  as  a  potential  means  for  fourth-order  accurate 
solutions  of  Navier-Stokes  equations  and  their  approximation  forms. 

To  improve  accuracy  and  simultaneously  achieve  superior  convergence, 

U.  Ghia,  K.  Ghia  and  Shin  (1981)  developed  a  multi-grid  (MG)  method  for  the 
solution  of  the  two-dimensional  Navier-Stokes  equations.  A  coupled  strongly 
implicit  multi-grid  (CSI-MG)  method  was  used  for  simultaneous  solution  of 
the  vorticity  and  stream  function  equations.  The  model  problem  used  was 
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that  of  the  flow  inside  a  shear-driven  cavity.  The  potential  of  the  method 
was  demonstrated  via  efficient  computation  of  solutions  for  Reynolds  number 
4 

up  to  10  using  a  very  fine  grid.  Because  of  the  appearance  of  one  or  more 
secondary  vortices  in  the  flow  field,  uniform  mesh  refinement  was  preferred 
to  the  use  of  one-dimensional  grid  clustering  coordinate  transformations. 
This  method  was  further  extended  by  K.  Ghia,  U.  Ghia  and  Shin  (1981)  for  use 
in  determining  asymptotic  three-dimensional  flow  in  curved  ducts  and  by 
U.  Ghia,  K.  Ghia  and  Ramamurti  (1983)  for  determining  the  developing  three- 
dimensional  flow  in  curved  ducts.  In  this  latter  parabolized  Navier-Stokes 
analysis,  the  MG  method  was  advanced  for  use  with  Neumann  boundary-value 
problems  in  clustered  curvilinear  coordinates.  This  comprised  an  important 
step  in  the  analysis  of  viscous  flows  using  the  velocity-pressure 
formulation  of  the  Navier-Stokes  equations.  With  successive  over-relaxation 
(SOR)  as  the  smoothing  operator  and  with  suitably  formulated  restriction  and 
coarse-grid-correction  operators,  a  4-grid  MG  procedure  enhanced  the 
efficiency  of  fine-grid  solutions  of  the  Neumann  problem  approximately  by  a 
factor  of  four.  U.  Ghia,  K.  Ghia  and  Ramamurti  (1983)  also  carefully 
examined  the  influence  of  the  smoothing  operator  by  employing  the  ADI  and  SI 
techniques  in  place  of  SOR.  For  the  curved  polar  duct  with  a  clustered 
grid,  Fig.  10  shows  the  convergence  history  of  the  MG-SI  method  for  the 
Neumann  Poisson  problem.  The  computational  advantage  of  the  MG  procedure 
with  increasing  refinement  of  the  finest  grid  is  very  obvious  from  this 
figure. 

As  stated  earlier  in  this  section,  the  MG-SI  method  was  applied  to  the 
shear-driven  cavity  problem.  The  MG-SI  method  developed  by  U.  Ghia,  K.  Ghia 
and  Shin  (1982)  had  about  fourteen  times  faster  convergence  rate  as  compared 
to  the  available  solvers,  including  that  of  Benjamin  and  Denny  (1979).  The 
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detailed  structure  of  the  various  secondary  eddies  at  Re  =  10  is  shown  in 
Fig.  11a;  these  have  become  the  benchmark  solutions  for  this  fairly  standard 
model  problem  as  they  were  obtained  using  a  fine  grid  of  (257*257)  points. 
Further,  K.  Ghia  and  U.  Ghia  (1984)  suggested,  via  the  results  presented  in 
Fig.  11b,  that  the  original  experiments  of  Pan  and  Acrivos  (1967)  were 
three-dimensional  in  nature.  Koseff  and  Street  (1983)  repeated  the 
experiment  with  a  spanwise  aspect  ratio  of  three  and  confirmed  that  the 
predictions  of  U.  Ghia,  K.  Ghia  and  Shin  (1982)  are  correct;  these  new 
experimental  results  are  also  shown  in  Fig.  11b. 

The  convergence  of  the  conventional  ADI  method  was  also  accelerated  by 
using  a  local  time-step  method  suggested  by  K.  Ghia  (1975).  Abdelhalim, 

U.  Ghia  and  K.  Ghia  (1983)  found  that  it  was  almost  essential  to  use  this 
procedure  to  obtain  converged  results  for  their  problem  of  flow  past  a  blunt 
plate  including  a  separation  bubble.  The  approach  can  be  viewed  as  a  means 
of  preconditioning  the  iteration  matrix  of  the  relaxation  scheme  and 
corresponds  to  a  more  uniform  Courant  number  throughout  the  flow  field. 

2.6  Fully-Impllcit  Numerical  Methods;  Semi-Direct  and  Direct  Methods 

The  development  of  semi-direct  and  direct  Poisson  solvers  was  also 
given  considerable  attention.  The  fast  direct  solvers  of  Hockney  (1970)  and 
Buneman  (1969)  were  widely  used  for  two-dimensional  flow  problems,  involving 
separable  Poisson  equations.  It  was  felt  that  these  methods  were  not  easily 
extendable  to  the  generalized  Poisson  equation  or  to  other  generalized 
elliptic  equations.  Therefore,  initially,  a  Cauchy-R iemann  solver  which 
leads  to  a  semi-direct  (SD)  method  of  solution  was  developed.  This  method 
followed  that  of  Martin  (1978)  such  that  the  form  of  the  generalized 


operator  in  non-orthogonal  coordinates  was  made  to  fit  the  form  of  the 


Cauchy-R iemann  operator  in  rectangular  coordinates.  The  results  of  this 
feasibility  study  were  given  by  Osswald,  K.  Ghia  and  U.  Ghia  (1980).  These 
results  showed  that  the  degree  of  non-orthogonality  and  grid  clustering 
strongly  influenced  the  solution  convergence  rate.  However,  it  remains  a 
viable  approach  for  efficient  solutions  of  the  Poisson  equation  in  a  variety 
of  orthogonal  and  non-orthogonal  coordinate  systems. 

For  solution  of  the  Poisson  equation  in  generalized  orthogonal 
coordinates,  Osswald  and  K.  Ghia  (1981)  developed  a  direct  block  Gaussian 
elimination  (BGE)  method.  Block-Gaussian  elimination  is  a  direct  extension 
of  the  Gaussian  elimination  procedure  to  block  matrices.  The  BGE  procedure 
provides  the  effective  inversion  of  an  [(NM)x(NM)]  matrix  through  the  actual 
inversion  of  a  predetermined  sequence  of  N(MxM)  sub-matrices.  The  BGE 
method  was  compared  with  the  SD  method  of  Martin  (1978)  for  the  general 
Poisson  problem  for  evaluation  of  accuracy  and  efficiency  and  was  found  to 
yield  a  direct  one-step  solution,  irregardless  of  the  degree  of  grid 
clustering,  with  considerably  increased  efficiency  as  compared  to  the  SD 
method.  A  comparison  of  various  existing  subroutines  for  solving  the  two- 
dimensional  Dirichlet  Poisson  problem  was  also  provided  and,  as  shown  in 
Table  1 ,  direct  BGE  provides  a  solver  with  efficiency  comparable  to  other 
available  direct  solvers,  but  with  the  added  advantage  that  it  is  applicable 
for  unsteady  flow  analysis  in  generalized  orthogonal  coordinates.  This  BGE 
solver  was  also  extended  by  Osswald,  K.  Ghia  and  U.  Ghia  (1983)  for  solution 
of  the  Neumann  Poisson  problem  encountered  for  the  pressure  field  in  an 
unsteady  pr imiti ve-var iable  Navier-Stokes  analysis.  The  extension  of  the 
BGE  method  to  the  solution  of  the  Poisson  equation  for  axisymmetric  flow  was 
also  provided  by  Osswald,  K.  Ghia  and  U.  Ghia  ( 1 98U ) ,  whereas  the  solution 
of  the  Poisson  equation  for  a  doubly-inf inite  region  was  provided  by 
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Of  the  various  CFD  analyses  developed,  some  were  of  direct  use  to  the 


technical  community.  Although  to  our  knowledge,  none  of  these  analyses  were 
used  in  the  development  of  any  specific  hardware,  they  are  being  used  in 
preliminary  design  studies  by  analysts  in  the  industry.  Some  of  these 
analyses  are  also  being  used  by  other  researchers  at  governmental 
laboratories  to  improve  their  analyses.  The  following  is  a  list  of  the  CFD 
analyses  and  the  organizations  using  them. 


ANALYSIS 

©  The  Parabolized  Navier-Stokes  Analysis 
for  Three-Dimensional  Internal  Flows 

©  The  Unsteady  Navier-Stokes  Analysis 
for  Internal  Flows 

©  The  Unsteady  Navier-Stokes  Analysis 
for  External  Flows 

©  Adaptive  Grid  Generation  Analysis 
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General  Electric  Co.,  Cincinnati,  OH; 
AVCO  Corp.,  Everett,  MA. 

General  Electric  Co.,  Cincinnati,  OH. 


McDonnell  Aircraft,  Co.,  St.  Louis,  MO; 
NASA-Langley  Research  Cntr.,  Hampton,  VA. 

Sverdrup  Technology,  Inc.,  Arnold 
Air  Force  Station,  TN. 


©  Navier-Stokes  Analysis  in  Primitive 
Variables 


NASA  Langley  Research  Center, 
Hampton,  VA. 
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ABLE  1.  COMPARISON  OF  SUBROUTINES  FOR  THE  2-D  DIRICHLET  POISSON  PROBLEM 
FOLLOWING  SCHUMANN  (1977). 


STREAMLINE  CONTOURS. 


ADIABATIC  WALL 


IG.  2a.  COMPARATIVE  STUDY  OF  ISOTROPIC  AND  ANISOTROPIC  TURBULENT  MODELS, 


INJECTION  P6(u=.995) 


b  (ii)  STEADY-STATE  SOLUTION  (AT  T=4l)  WITH  IMPROVED  MESH, 


FIG.  3, 


STREAM-FUNCTION  CONTOURS  FOR  CONSTRICTED  CHANNEL  FLOW, 
Re  =  1000,  Ai|i  =  0.002  WITHIN  SEPARATION  BUBBLE; 

A* =  0.1  OTHERWISE. 
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COMPARISON  OF  MEASURED  AND  COMPUTED  RECIRCULATiON  LENGTHS  FOR  BLUNT 
PLATES  AND  LONGITUDINAL  CYLINDERS. 


DIFFERENTIAL  EQUATION;  u  -  (u-U) Re  u  =  Re  u  for  -5  <  x  <  5  and 


SOLUTION  OF  1-D  VISCOUS  BURGERS  EQUATION,  WITH  RE  RANGING  FROM  ]  TO  10,000 


29  49  69  89  109 

FINE-GRID  ITERATIONS 


CONVERGENCE  HISTORY  OP  MULTIGRID  METHOD  FOR  NEUMANN 
PRESSURE  PROBLEM  IN  DEVELOPING  FLOW  THROUGH  CURVED 
POLAR  DUCT  WITH  GRID  CLUSTERING, 


